Distribution of the number of significant effect sizes

A while back, I posted the outline of a problem about the number of significant effect size estimates in a study that reports multiple outcomes. This problem interests me because it connects to the issue of selective reporting of study results, which creates problems for meta-analysis. Here, I’ll re-state the problem in slightly more general terms and then make some notes about what’s going on.

Consider a study that assesses some effect size across $m$ different outcomes. (We’ll be thinking about one study at a time here, so no need to index the study as we would in a meta-analysis problem.) Let $T_i$ denote the effect size estimate for outcome $i$, let $V_i$ denote the sampling variance of the effect size estimate for outcome $i$, and let ${\theta }_i$ denote the true effect size parameter for corresponding to outcome $i$. Assume that the study outcomes ${\left[T_i\right]}_{i=1}^m$ follow a correlated-and-hierarchical effects model, in which $$T_i=\mu +u+v_i+e_i,$$ where the study-level error $u\sim N(0,{\tau }^2)$, the effect-specific error $v_i\stackrel{iid}{\sim }N(0,{\omega }^2)$, and the vector of sampling errors ${\left[e_i\right]}_{i=1}^m$ is multivariate normal with mean $\mathbf{0}$, known variances $\text{Var}(e_i)={\sigma }^2$, and compound symmetric correlation structure $\text{cor}(e_h,e_i)=\rho$.

Define $A_i$ as an indicator that is equal to one if $T_i$ is statistically significant at level $\alpha$ based on a one-sided test, and otherwise equal to zero. (Equivalently, let $A_i$ be equal to one if the effect is statistically significant at level $2\alpha$ and in the theoretically expected direction.) Formally, $$A_i=I(\frac{T_i}{\sigma }>q_{\alpha })$$ where $q_{\alpha }={\mathrm{\Phi }}^{-1}(1-\alpha )$ is the critical value from a standard normal distribution (e.g., $q_{.05}=1.645$, $q_{.025}=1.96$). Let $N_A=\sum _{i=1}^m{A}_i$ denote the total number of statistically significant effect sizes in the study. The question is: what is the distribution of $N_A$.

Compound symmetry to the rescue

As I noted in the previous post, this set-up means that the effect size estimates have a compound symmetric distribution. We can make this a bit more explicit by writing the sampling errors in terms of the sum of a component that’s common acrosss outcomes and a component that’s specific to each outcome. Thus, let $e_i=f+g_i$, where $f\sim N(0,\rho {\sigma }^2)$ and $g_i\stackrel{iid}{\sim }N(0,(1-\rho ){\sigma }^2)$. Let me also define $\zeta =\mu +u+f$ as the conditional mean of the effects. It then follows that the effect size estimates are conditionally independent, given the common components: $$(T_i|\zeta )\stackrel{iid}{\sim }N(\zeta ,{\omega }^2+(1-\rho ){\sigma }^2)$$ Furthermore, the conditional probability of a significant effect is $$\text{Pr}(A_i=1|\zeta )=\mathrm{\Phi }\left(\frac{\zeta -q_{\alpha }\sigma }{\sqrt{{\omega }^2+(1-\rho ){\sigma }^2}}\right)$$ and $A_1,...,A_m$ are mutually independent, conditional on $\zeta$. Therefore, the conditional distribution of $N_A$ is binomial, $$(N_A|\zeta )\sim Bin(m,\pi )$$ where $$\pi =\mathrm{\Phi }\left(\frac{\zeta -q_{\alpha }\sigma }{\sqrt{{\omega }^2+(1-\rho ){\sigma }^2}}\right).$$ What about the unconditional distribution?

To get rid of the $\zeta$, we need to integrate over its distribution, which leads to $$\text{Pr}(N_A=a)=\text{E}\left[\text{Pr}(N_A|\zeta )\right]=\int f_{N_A}(a|\zeta ,\omega ,\sigma ,\rho ,m)\times f_{\zeta }(\zeta |\mu ,\tau ,\sigma ,\rho )d\zeta ,$$ where $f_{N_A}(a|\zeta ,\omega ,\sigma ,\rho )$ is a binomial density with size $m$ and probability $\pi =\pi (\zeta ,\omega ,\sigma ,\rho )$ and $f_{\zeta }(\zeta |\mu ,\tau ,\sigma ,\rho )$ is a normal density with mean $\mu$ and variance ${\tau }^2+\rho {\sigma }^2$.

This distribution is what you might call a binomial-normal convolution or a random-intercept probit model (where the random intercept is $\zeta$). As far as I know, the distribution cannot be evaluated analytically but instead must be calculated using some sort of numerical integration routine.

Just the moments, please

If all we care about is the expectation of $N_A$, we don’t need to bother with all the conditioning business and can just look at the marginal distribution of the effect size estimates taken individually. Marginally, $T_i$ is normally distributed with mean $\mu$ and variance ${\tau }^2+{\omega }^2+{\sigma }^2$, so $\text{Pr}(A_i=1)=\psi$, where $$\psi =\mathrm{\Phi }\left(\frac{\mu -q_{\alpha }\sigma }{\sqrt{{\tau }^2+{\omega }^2+{\sigma }^2}}\right).$$ By the linearity of expectations, $$\text{E}(N_A)=\sum _{i=1}^m\text{E}(A_i)=m\psi .$$

We can also get an approximation for the variance of $N_A$ by working with its conditional distribution above. By the rule of variance decomposition, $$\begin{array}{rl}\text{Var}(N_A)& =\text{E}\left[\text{Var}(N_A|\zeta )\right]+\text{Var}\left[\text{E}(N_A|\zeta )\right]\\ & =m\times \text{E}[\pi (1-\pi )]+m^2\times \text{Var}\left[\pi \right]\\ & =m\times \text{E}\left[\pi \right](1-\text{E}\left[\pi \right])+m(m-1)\times \text{Var}\left[\pi \right],\end{array}$$ where $\pi$ is, as defined above, a function of $\zeta$ and thus a random variable. Now, $\text{E}(\pi )=\psi$ and we can get something close to $\text{Var}(\pi )$ using a first-order approximation: $$\text{Var}\left(\pi \right)\approx {\left({\frac{\delta \pi }{\delta \zeta }|}_{\zeta =\mu }\right)}^2\times \text{Var}\left(\zeta \right)={\left[\varphi \left(\frac{\mu -q_{\alpha }\sigma }{\sqrt{{\omega }^2+(1-\rho ){\sigma }^2}}\right)\right]}^2\times \frac{{\tau }^2+\rho {\sigma }^2}{{\omega }^2+(1-\rho ){\sigma }^2}.$$ Thus, $$\begin{array}{r}\text{Var}(N_A)\approx m\times \psi (1-\psi )+m(m-1)\times {\left[\varphi \left(\frac{\mu -q_{\alpha }\sigma }{\sqrt{{\omega }^2+(1-\rho ){\sigma }^2}}\right)\right]}^2\times \frac{{\tau }^2+\rho {\sigma }^2}{{\omega }^2+(1-\rho ){\sigma }^2}.\end{array}$$ If the amount of common variation is small, so ${\tau }^2$ is near zero and $\rho$ is near zero, then the contribution of the second term will be small, and $N_A$ will act more or less like a binomial random variable with size $m$ and probability $\psi$. On the other hand, if the amount of independent variation in the effect sizes is small, so ${\omega }^2$ is near zero and $\rho$ is near 1, then the term on the right will approach $m(m-1)\psi (1-\psi )$ and $\text{Var}\left(N_A\right)$ will approach $m^2\psi (1-\psi )$, or the variance of $m$ times a single Bernoulli variate. So you could say that $N_A$ has anywhere between $1$ and $m$ variate’s worth of information in it, depending on the degree of correlation between the effect size estimates.

Interactive distribution

Here is an interactive graph of the probability mass function of $N_A$, with probability points calculated using Gaussian quadrature. Below the graph, I also report $\psi$, the exact mean and variance of $N_A$, and the first-order approximation to the variance (denoted $V_{approx}). When $\tau >0$ and $\rho >0$, the approximate variance is not all that accurate because the first-order approximation to $\text{Var}(\pi )$ isn’t that good.

math = require("mathjs")
norm = import('https://unpkg.com/norm-dist@3.1.0/index.js?module')

quad_points = JSON.parse(all_quad_points).at(qp - 1)

sigma = 2 / math.sqrt(ESS)

zeta_sd = math.sqrt(tau**2 + rho * sigma**2)
ID_sd = math.sqrt(omega**2 + (1 - rho) * sigma**2)

crit = norm.icdf(1 - alpha)

binomial_coefs = Array(m+1).fill(null).map((x,index) => {
  return math.combinations(m, index);
})

probs = quad_points.map(zeta => {
  let Z = (zeta[0] * zeta_sd + mu - crit * sigma) / ID_sd;
  return [norm.cdf(Z), zeta[1]];
})

p_binom_norm = binomial_coefs.map((coef, a) => {
  let p = probs.map((x) => {
    return (x[0]**a) * ((1 - x[0])**(m - a)) * x[1];
  });
  return coef * math.sum(p);
})

math = Object {isNumber: ƒ(e), isComplex: ƒ(e), isBigNumber: ƒ(e), isBigInt: ƒ(e), isFraction: ƒ(e), isUnit: ƒ(e), isString: ƒ(e), isArray: ƒ(), isMatrix: ƒ(e), isCollection: ƒ(e), isDenseMatrix: ƒ(e), isSparseMatrix: ƒ(e), isRange: ƒ(e), isIndex: ƒ(e), isBoolean: ƒ(e), isResultSet: ƒ(e), isHelp: ƒ(e), isFunction: ƒ(e), isDate: ƒ(e), isRegExp: ƒ(e), …}

norm = Module {Z: ƒ(…), cdf: ƒ(z), icdf: ƒ(…), intE: ƒ(a, b), pdf: ƒ(z), Symbol(Symbol.toStringTag): "Module"}

quad_points = Array(21) [Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), …]

sigma = 0.22360679774997896

zeta_sd = 0.19999999999999998

ID_sd = 0.17320508075688773

crit = 1.959963986120195

binomial_coefs = Array(7) [1, 6, 15, 20, 15, 6, 1]

probs = Array(21) [Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), Array(2), …]

p_binom_norm = Array(7) [0.3632983110161876, 0.18642239532408006, 0.13279056012006735, 0.10429406708513911, 0.0848789783106351, 0.0702147931381166, 0.058003374919324134]

psi = norm.cdf((mu - crit * sigma) / math.sqrt(tau**2 + omega**2 + sigma**2))
psi_print = psi.toFixed(3)

E_NA = m * psi 
E_NA_print = E_NA.toFixed(3)

dpi_dzeta = norm.pdf((mu - crit * sigma) / ID_sd)
V_pi_approx = (dpi_dzeta * zeta_sd / ID_sd)**2
V_approx = m * psi * (1 - psi) + m * (m - 1) * V_pi_approx
V_approx_print = V_approx.toFixed(3)

V_NA = {
  let V_NA = 0 - E_NA**2;
  for (let i = 0; i <= m; i++) {
    V_NA += i**2 * p_binom_norm[i];
  }
  return V_NA;
}
V_NA_print = V_NA.toFixed(3)

psi = 0.3006337901985208

psi_print = "0.301"

E_NA = 1.803802741191125

E_NA_print = "1.804"

dpi_dzeta = 0.290096539705762

V_pi_approx = 0.11220800313234232

V_approx = 4.627758780306626

V_approx_print = "4.628"

V_NA = 3.60408188896073

V_NA_print = "3.604"

Distribution of $N_A$

Plot.plot({
  x: {
    label: "Number of significant effect sizes"
  },
  y: {
    domain: [0, 1],
    label: "Probability"
  },
  marks: [
    Plot.ruleY(0),
    Plot.barY(p_binom_norm, {
      fill: "steelblue"
    }),
  ]
})

Moments of $N_A$

tex`
\begin{aligned}
\psi &= ${psi_print} \\
\mathbb{E}\left(N_A\right) &= ${E_NA_print} \\
\mathbb{V}\left(N_A\right) &= ${V_NA_print} &V_{approx} &= ${V_approx_print}
\end{aligned}
`

$\begin{aligned} \psi &= 0.301 \\ \mathbb{E}\left(N_A\right) &= 1.804 \\ \mathbb{V}\left(N_A\right) &= 3.604 &V_{approx} &= 4.628 \end{aligned}$

viewof m = Inputs.range(
  [1, 30], 
  {value: 6, step: 1, label: "Number of effect sizes (m):"}
)

viewof ESS = Inputs.range(
  [4, 300], 
  {value: 80, step: 1, label: "Effective sample size:"}
)

viewof mu = Inputs.range(
  [-2, 2], 
  {value: 0.3, step: 0.01, label: "Average effect size (mu):"}
)

viewof tau = Inputs.range(
  [0, 1], 
  {value: 0.1, step: 0.01, label: "Between-study SD (tau):"}
)

viewof omega = Inputs.range(
  [0, 1], 
  {value: 0.1, step: 0.01, label: "Within-study SD (omega):"}
)

viewof rho = Inputs.range(
  [0, 1], 
  {value: 0.6, step: 0.01, label: "Sampling error correlation (rho):"}
)

viewof alpha = Inputs.range(
  [0.005, 0.995], 
  {value: 0.025, step: .005, label: "One-sided significance threshold (alpha):"}
)

viewof qp = Inputs.range(
  [1, 30], 
  {value: 21, step: 1, label: "Number of quadrature points:"}
)

Number of effect sizes (m):

m = 6

Effective sample size:

ESS = 80

Average effect size (mu):

mu = 0.3

Between-study SD (tau):

tau = 0.1

Within-study SD (omega):

omega = 0.1

Sampling error correlation (rho):

rho = 0.6

One-sided significance threshold (alpha):

alpha = 0.025

Number of quadrature points:

qp = 21